Early signatures of regime shifts in gene expression dynamics 
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Recently, a large number of studies have been carried out on the early signatures of sudden 
regime shifts in systems as diverse as ecosystems, financial markets, population biology and complex 
diseases. Signatures of regime shifts in gene expression dynamics are less systematically investigated. 
In this paper, we consider sudden regime shifts in the gene expression dynamics described by a fold- 
bifurcation model involving bistability and hysteresis. We consider two alternative models, Models 1 
and 2, of competence development in the bacterial population B. subtilis and determine some early 
signatures of the regime shifts between competence and noncompetence. We use both deterministic 
and stochastic formalisms for the purpose of our study. The early signatures studied include the 
critical slowing down as a transition point is approached, rising variance and the lag-1 autocorrelation 
function, skewness and a ratio of two mean first passage times. Some of the signatures could 
provide the experimental basis for distinguishing between bistability and excitability as the correct 
mechanism for the development of competence. 

PACS numbers: 87.10.Mn, 87.16. Yc, 87.17.Aa, 87.18.Cf, 87.18.Tt 



I. INTRODUCTION 

Complex dynamical systems, ranging from ecosystems 
and the climate to gene regulatory networks and financial 
markets, are known to exhibit abrupt shifts from one dy- 
namical regime to a contrasting one at critical parameter 
values Examples of sudden regime shifts include 

the collapse of vegetation when semi-arid conditions pre- 
vail, the transition from a clear lake to a turbid one, sud- 
den changes in fish/ wildlife populations [1, distinct 
changes in the climate and patterns of oceanic circula- 
tion [l|, Q, financial markets undergoing global crashes 
1|, spontaneous systemic failures such as asthma attacks 
7 1, or epileptic seizures [|[ etc. In a gene regulatory net- 
work, a sudden transition may occur from one stable gene 
expression state to another at a critical parameter value 
[9l4l2l|. The induction of the lac operon in E. coli results 
in the synthesis of the protein /3-galactosidase required 
for breaking up sugar molecules and releasing energy to 
the cell. There is now experimental evidence Q that 
an abrupt transition from the uninduced ((3-galactosidase 
level low) state to the induced (fi-galactosidase level high) 
state of the lac operon occurs at a critical value of an in- 
ducer concentration. 

The regime/state shifts in the examples mentioned 
above are mostly a consequence of the fold-bifurcation 
(or the fold-catastrophe) , well-characterized in dynamical 
systems theory [HUH, El- Figure [T] illustrates a specific 
type of the fold-bifurcation based on bistability and hys- 
teresis, which provides a physical understanding of the 
features associated with sudden regime shifts. The plot 
represents the steady states of a dynamical system versus 
a specific parameter. The state of the dynamical system 
at time t is defined by the magnitudes of the relevant 
variables at that time. In the steady state, the net rate of 
change in the magnitude of a variable is zero so that there 
is no further time evolution. The solid lines in Figure [T] 
denote stable steady states separated by a branch (dotted 
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line) of unstable steady states. The steady state curve is 
folded backwards giving rise to bistability, i.e., the exis- 
tence of two stable steady states in the shaded region. 
The parameter values u\ and ui represent the bifurca- 
tion points at which the abrupt regime shifts from one 
stable steady state to another occurs. On crossing the bi- 
furcation point, the system loses bistability and becomes 
monostable. The transition from one branch of stable 
steady states to the other is not reversible but exhibits 
hysteresis. This implies that if a transition takes place at 
a bifurcation point, say the upper one, the reverse tran- 
sition from the upper to the lower branch occurs at the 
lower bifurcation point and not at the upper bifurcation 
point itself. Bistability owes its origin to the presence of 
one or more positive feedback loops governing the system 
dynamics with the added condition that the dynamics be 
sufficiently nonlinear [Tol - iT^ | . Each stable steady state 
has its own basin of attraction in the state space, i.e., the 
space of all states. In the case of more than one stable 
steady state, the location of the initial state of the sys- 
tem in a basin of attraction determines the steady state 
attained by the system in the course of time. 

The stability of a steady state indicates that the sys- 
tem returns to the state on being weakly perturbed from 
it. This is shown by arrow directions in Figure [TJ If 
the perturbation is sufficiently strong so that a transition 
takes place from one basin of attraction to the other, a 
switch occurs between the stable steady states even be- 
fore the bifurcation point is reached. Closer the system 
is to the bifurcation point, the lesser is the magnitude of 
the perturbation needed to bring about the switch. In 
the example of Figure [TJ the branch of unstable steady 
states constitutes the border between the two basins of 
attraction. The distance of a stable steady state from 
the corresponding unstable steady state is a measure of 
the resilience (robustness against perturbation) of the sys- 
tem. The resilience is gradually weakened as the system 
approaches a bifurcation point. The dynamics of natu- 
ral systems, in general, have a stochastic component due 
to the presence of random external influences and the in- 
herently probabilistic nature of the processes involved in 
the dynamics. Consider the dynamics of gene expression 
in a gene regulatory network. Gene expression consists 
of two major steps: transcription and translation during 



Bifurcation Parameter 

Figure 1. A generic steady state versus bifurcation parameter 
diagram. The shaded region represents the region of Insta- 
bility separating two regions of monostability. The solid lines 
correspond to stable steady states and the dashed line repre- 
sents unstable steady states. The points u\ and U2 are the 
lower and upper bifurcation points respectively. The arrows 
indicate the time evolution of a weakly perturbed system with 
the system moving towards stable steady states and moving 
away from unstable steady states 



which messenger RNA (mRNA) and protein molecules re- 
spectively are synthesized. The biochemical events (e.g., 
the binding of a RNA polymerase at the promoter re- 
gion of the DNA to initiate transcription) constituting 
gene expression are inherently stochastic in character re- 
sulting in fluctuations (noise) around mean mRNA and 
protein levels [3, UH . Extrinsic influences like variability 
in the number of regulatory molecules also contribute to 
the noise. Instead of a single steady state protein level, as 
in the case of deterministic time evolution, one now has a 
steady state probability distribution in the protein levels 
reflecting the stochastic nature of the time evolution. In 
the presence of low/moderate amounts of noise, the phys- 
ical picture underlying sudden regime shifts still remains 
valid. As in the case of applied perturbations, fluctuations 
of sufficiently strong magnitude can bring about regime 
shifts before the bifurcation point is crossed. A number of 
early signatures of regime shifts have been proposed so far 
0} EL EH, E3 m the scenario of the fold-bifurcation. These 
are the critical slowing down (CSD) and its associated ef- 
fects, namely, rising variances and the lag-1 autocorrela- 
tion function as the critical transition point is approached, 
increased skewness in the steady state probability distri- 
bution and the presence of flickering transitions between 
the stable steady states. The utility of such signals in 
the cases of ecological and financial systems [l| and com- 
plex deteriorating diseases [l8| is significant, specially, in 
developing appropriate risk-aversion/management strate- 
gies. Also, in the absence of physical models capturing the 
essential dynamics, one can obtain quantitative measures 
of the early signatures by analyzing time-series data. 

In this paper, we compute quantities providing early 
signatures of abrupt regime shifts in gene expression dy- 
namics. The dynamics undergo the fold-bifurcation re- 
sulting in binary gene expression, i.e., the existence of two 
stable expression states. We specifically focus on the phe- 
nomenon of competence development in the soil bacteria 
B. subtilis as a manifestation of binary gene expression. 




(b) 



Figure 2. Two models, Model 1 (a) and Model 2 (b) of com- 
petence development, (a) The comK gene expresses ComK 
proteins (represented by A) which form dimers. A dimer binds 
the promoter region of the gene and activates the initiation of 
transcription. The dimer can also unbind from the promoter 
region and dissociate into free monomers, (b) The autoregu- 
latory positive feedback loop mediated by ComK proteins is 
present. In addition there is a negative feedback loop in which 
ComK inhibits the expression of the comS gene and the ComS 
proteins repress the activity of the MecA-ClpP-ClpC Complex 
which targets ComK proteins for degradation. 



Microorganisms like bacteria have to cope with a num- 
ber of stresses during their lifetime. The bacteria adopt 
a number of strategies to enhance their chances of sur- 
vival under changed circumstances [H, EH ■ One of these 
is the generation of phenotypic heterogeneity (no genetic 
change) in the bacterial population so that a fraction of 
the population, even if small, may survive under stressful 
conditions. In the B. subtilis population, only a fraction 
of the population, in which the level of a key transcription 
factor ComK is high, develops competence [HI, [2(| . The 
protein acts as a master regulator activating the tran- 
scription of several genes including those essential for the 
uptake of foreign DNA from the environment. The incor- 
poration of the new DNA into the bacterial genome con- 
fers favorable traits on the bacterial subpopulation with 
high ComK level, enabling the subpopulation to adapt to 
stress. The ComK activity resulting in the development 
of competence is in turn controlled by a host of other pro- 
teins. The core module of the complex regulatory network 
consists of an autoregulatory positive feedback loop in 
which the ComK proteins promote their own production. 
The positive feedback gives rise to binary gene expression 
in the cell population, i.e., two distinct subpopulations 
with low and high ComK levels respectively. Two inde- 
pendent experiments [2ll.[22j| provide confirmation that an 
autoregulatory positive feedback loop of ComK produc- 
tion is by itself sufficient to establish binary gene expres- 
sion in a bacterial population. Considering deterministic 
time evolution, the steady state ComK level versus an ap- 
propriate gene expression parameter exhibits a hysteresis 
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curve, resulting from the fold-bifurcation, similar to the 
one shown in Figure [Tl [23l| . In this scenario, if the cells in 
a population are prepared to be in the same initial state, 
all the cells evolve to the same final state giving rise to 
a homogeneous cell population. The generation of het- 
erogeneity in the form of two distinct subpopulations is 
brought about by fluctuation-driven transitions between 
the low and high ComK expression states, the fluctua- 
tions being associated with the ComK levels. This gives 
rise to the experimentally observed [HI, [2(| bimodal dis- 
tribution in the ComK levels, i.e., a distribution with two 
prominent peaks. There is now experimental evidence 
[24| that a lower fraction of the B. subtilis population de- 
velops competence with reduced noise in the low ComK 
level. An alternative physical mechanism underlying com- 
petence development has been proposed by Siiel et. al. 
[2 51 ] in terms of excitability in the dynamics of the ge- 
netic circuit. The excitable core module consists of both 
positive and negative feedback loops which bring about a 
transient activation to the high ComK state. In this case, 
there is only one stable steady state (low ComK level) and 
two unstable steady states the lower of which, in terms of 
the ComK level sets a threshold for switching j25|. Fluc- 
tuations in the low ComK level activate the switch to ex- 
pression states in the neighborhood of the state with high 
ComK level which constitutes an unstable steady state. 
The transient activation is followed by an ultimate return 
to the stable low ComK state. At any instant of time, 
the population divides into two subpopulations with low 
and high ComK levels respectively, signifying a different 
origin for the bimodal distribution. Quantitative time 
lapse fluorescence microscopy provides experimental evi- 
dence [25| of the probabilistic and transient differentiation 
into competence. In this paper, we study sudden regime 
shifts in the gene expression dynamics leading to com- 
petence development in B. subtilis We provide a quanti- 
tative characterization of the early signatures of regime 
shifts, both noise-driven and occurring at the bifurcation 
points, in the fold-bifurcation model. In Section UH we 
describe two models of competence development, Models 
1 and 2, exhibiting bistable and excitable dynamics. A 
stochastic component is added to only Model 1 which is 
a one- variable model. Model 2 involves two variables and 
is treated dctcrministically. The dynamics of the stochas- 
tic version of Model 1 are described by the Langevin and 
Fokkcr-Planck equations. We further define the quanti- 
tative measures of the early signatures of sudden regime 
shifts. In Section ITO1 we present the results of our com- 
putations on these early signatures. Section IIVI contains 
concluding remarks. 



II. MODELS AND EARLY SIGNATURES OF 
REGIME SHIFTS 

We first consider a simple model (Model 1) in which the 
protein product of a single gene autoactivates its own pro- 
duction via a positive feedback loop. The genetic circuit 
is shown in Figure [U^a) and, as mentioned earlier, con- 
stitutes a core module of the complex network resulting 
in competence development. The proteins synthesized by 
the comK gene form dimers. The dimer molecules bind 
at the promoter region of the DNA and activate gene ex- 



pression. The gene also synthesizes proteins at a basal 
level. A kinetic scheme of the model is as follows 1231: 



P 2 +G 



GP 2 ^ G* 



G^P 



G* 
P + P 



P 

P 2 



pH 



(1) 



The gene can be in two possible states: inactive (G) and 
active (G*). Proteins are synthesized with rate constant 
Ji(Jo) in the state G*{G) with Jo <C Ji- The synthesized 
proteins dimerize with K e as the equilibrium dissociation 
constant. The protein dimer P 2 binds to the gene in state 
G and forms the complex GP 2 with k\ and k 2 being the 
rate constants for binding and unbinding. The complex 
GP 2 in turn is activated to the state G* , the rate constants 
k a and kd being the activation and deactivation rate con- 
stants. The synthesized proteins are degraded with a rate 
constant, k p , 4> denoting the degradation product. As 
shown by Karmakar and Bose [23j, the kinetic scheme 
displayed in Equation ([I]) can be mapped onto a simpler 
scheme 
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The effective activation and deactivation rate constants 
k' a (x) and k' d are: 

In Equation ([3]), x denotes the protein concentration and 



k s = J j^K e . In the simplified scheme, the rate of change 
of the protein concentrations: 



dx 



Jokd 



(4) 



The steady state condition 4| = yields three solutions 
in a specific parameter regime corresponding to two stable 
steady states separated by an unstable steady state. The 
model dynamics exhibit the fold-bifurcation of the type 
shown in Figure [1] The rate constants Jo, Ji, and k a serve 
as the bifurcation parameters. Figure H fb) shows the ge- 
netic circuit (Model 2) , proposed in Ref. [25[ , as governing 
the dynamics of competence development via excitability. 
The circuit contains the autorcgulatory positive feedback 
loop of Figure HJa). In addition there is a negative feed- 
back loop in which the ComK protein represses the ex- 
pression of the comS gene whereas the ComS protein in- 
hibits the degradation of ComK by the MecA-ClpP-ClpC 
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Figure 3. The k a — kd phase diagram of Model 1 showing a 
region of bistability separating two regions of monost ability. 
The stability landscape U(x) versus x for three values of the 
bifurcation parameter is also shown. As one moves toward the 
lower bifurcation point k a = 1.344, indicated by an arrow in 
the region of bistability, the basin of attraction corresponding 
to the high stable expression state becomes flatter. At the 
bifurcation point itself, the return time Tr, associated with 
the CSD, becomes zero. 



Figure 4. The steady state probability distribution P a t(x) 
of protein levels in the region of bistability. The maxima of 
P 3 t{x) at xi and £3 represent low and high expression states 
respectively. The points X2 and X4 denote the lower and up- 
per cut-off points of the probability distribution P"(x) corre- 
sponding to the high expression state. The inset shows the 
stochastic potential 4>f(x) versus x. 



complex. The repression of comS by ComK is however, 
not well established under wild-type expression conditions 
[20I ] . The dynamics of the model are described in terms 
of the differential equations HU : 



dK 



K 



dt ak k% + K n 



dS_ 
~dt 



S 



S 



l + (K/k!)P 1 + K + S 



(5) 



(0) 



The variables K and S denote the concentrations of the 
ComK and ComS proteins respectively, ak and bk rep- 
resent the basal and fully activated rates of ComK syn- 
thesis and ko is the ComK concentration needed for 50% 
activation. The Hill coefficients n and p are indicative 
of the cooperativities of ComK autoactivation and ComS 
repression respectively. The expression rate of ComS has 
maximal value b s and is half- maximal at K = k\. The 
non-linear degradation terms are a consequence of the 
MecA complex-mediated degradation with a competitive 
inhibition of the degradation by ComS. 

A simple stochastic version of the Model 1 has been 
studied in Ref. 23]. In this model, the only stochasticity 
considered is associated with the random transitions of 
the gene between the inactive (G) and active (G*) states. 
The rest of the processes undergo deterministic time evo- 
lution. The simple model yields bimodal protein distribu- 
tions in the steady state in certain parameter regions. In 
this paper, the stochastic dynamics of the model are inves- 
tigated using the formulations based on the Langevin and 
Fokkcr-Planck (FP) equations. The steady state analysis 
of a bistable system in these formalisms is described in 
Refs. [26l - [3Cl |. The one- variable Langevin equation (LE) 
including both multiplicative and additive noise terms is 



given by 



(7) 



where e(t) (multiplicative noise) and T(t) (additive noise) 
are Gaussian white noises with mean zero and correla- 
tions: 

(e(t)e(tO) = 2d 1 S(t - t') 
(T(t)T(t')) =2d 2 S(t-t') 



(e(t)r(f)) = (T(t)e(t')) = 2X d ^d 2 5{t - t') (8) 

In Equation ©, d\ and d 2 denote the strengths of the 
noises e{t) and T(t) respectively and Xd is the degree of 
correlation between them. The first term in Equation 
([7]) describes the deterministic dynamics. In the case of 
Model 1, f(x) is given by the right hand side expression 
in Equation (U). The FP equation is a rate equation for 
the probability distribution P(x, t), obtained from the LE 
as (27L l35|: 

dPlx t) d d 2 

= - — [A(x)P(x,t)} + —[B(x)P(x 7 t)} (9) 

where 



A(x) = f(x) + d l9l (x)g[{ X ) + Xd^hg'^x) (10) 



and 



B(x)^d 1 [g 1 (x)} 2 +2X d ^hd 2gi {x)+d 2 (11) 
The steady state probability distribution (SSPD) is given 
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Figure 5. For Model 1, the variation of the variance a , the 
third moment M3, the coefficient of variation C„ and the mod- 
ulus I — y I where 7 is the skewness of the normalized steady state 
probability distribution P"(x) (Equation ((26j)) as the bifurca- 
tion parameter k a is decreased towards the lower bifurcation 
point. The solid lines correspond to the case when only addi- 
tive noise is present. The dotted lines are obtained when both 
additive and multiplicative noise terms are present in the LE 
(Equation (0). The strengths of the multiplicative and addi- 
tive noise are di = 0.25 and d,2 = 0.25 respectively. The other 



parameter values are kd 
and J 1 — 1. 
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where TV is the normalization constant. Equation (|12D 
can be rewritten in the form 



P st (x) = /Ve-<M*) 



(13) 



with 
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: f(y)dy 

di [gi (y)} 2 + 2A d V 'd\d 2 gi{y) 
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4>f{x) defines the "stochastic potential" of the dynamics. 
Once the steady state probability distribution is deter- 
mined, quantities like the variance and skewness, to be 
defined below, which provide early signatures of regime 
shifts can be determined. 

We next define the quantitative measures of the early 
signatures of regime shifts. The first signature is that of 
the CSD which is a distinguishing feature of the dynam- 
ics close to a bifurcation point [la, [l?], [H| . This involves 
a progressively larger relaxation time, as the bifurcation 



point is approached, to an attractor of the dynamics (say, 
a stable steady state) after being weakly perturbed from 
it. The physical origin of this phenomenon can be un- 
derstood in terms of the stability landscape U{x). The 
rate equation governing the dynamics of Model 1 (Equa- 
tion (g])) can be written as ^ = f(x) = Figure 
[3ja) shows the phase diagram of the model in the k a — k d 
plane with a region of bistability separating two regions 
of monostability. The other parameter values are J\ = 1, 
Jo = 0.01, k p = 0.03 in appropriate units. Figure \Mh) 
shows the three stability landscapes at three successively 
decreasing values of k a as one approaches the bifurca- 
tion point in the direction of the arrow from within the 
region of bistability. At point 1 (k a = 1.5), the states 
corresponding to the minima of U (x) represent the stable 
steady states. The location of the "ball" represents the 
state of the dynamical system. As one progresses from 
point 1 to point 2 (k a = 1.35), the steady state with 
high value of x becomes less stable. The associated basin 
of attraction becomes flatter with reduced size so that 
it takes a longer time for a perturbed state (ball shifted 
from the minimum position) to relax back to the stable 
steady state. The relaxation time diverges at the bifur- 
cation point where the stable steady state is on the verge 
of losing stability (point 3, k a = 1.344). A weak pertur- 
bation pushes the ball to the minimum with low value 
of x, i.e., the system does not relax back to the high x 
state. One can define a return time Tr which provides 
a measure of the time taken by the dynamical system to 
regain a stable steady state after being weakly perturbed 
from it. Let x s t be the stable steady state value of x and 
77(f) = x(t) — x st be the small deviation from the steady 
state value under weak perturbation. The deterministic 
rate equation is given by ^ = f(x). On Taylor expanding 
f(x) around x = x s t and keeping terms upto the order of 
77(f), one obtains 



dv . M . df(x) 
— = Xr){t), A = — — 
dt dx 



The solution of the Equation (|T5|) is 



7?(f) = 7,(0). 



(15) 



(16) 



where 77(0) is the initial value of 77(f) at time f = 0. 
The sign of the parameter A indicates the stability of the 
steady state. If A is < (> 0), the steady state is stable 
(unstable). We designate A as the stability parameter. It 
is a well-known result from dynamical systems theory that 
at a bifurcation point, the stability parameter A, associ- 
ated with the steady state losing stability, becomes zero 
EH, EE El- In the case of Model 1, one can check 
that A = at the two bifurcation points. The return time 
= m thus diverges at a bifurcation point. If the dy- 
namical system is described by more than one variable, 
the stability of a steady state is determined by the eigen- 
values of the Jacobian matrix governing the dynamics of 
the perturbed system [l3|. Let X max be the real part of 
the dominant eigenvalue of the Jacobian matrix (for sta- 
bility, the real parts of all the A's should be negative). The 
dominant eigenvalue is the one with the largest real part. 
The return time Tr is given by Tr = j^— — r. Examples 
of the experimental observations of the CSD include the 
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Figure 6. For Model 1, the plot of | 7 |, where 7 is the skew- 
ness of the normalized steady state probability distribution, 
versus the strength of additive noise d-2 for the parameter val- 
ues k a = 1.5, kd = 1.5, k p = 0.03, Jo = 0.01, Ji = 1, and 
di = (no multiplicative noise). The value k a = 1.5 falls in 
the region of bistability. 



transition from the G2 growth phase to the mitotic phase 
of the eukaryotic cell division cycle [3(|, a direct obser- 
vation of the CSD in a laboratory population of budding 
yeast before population collapse occurs at a critical ex- 
perimental condition |32j and the demonstration of the 
CSD in a population of cyanobacteria [33l ]. 

In the presence of noise, one can reexpress Equation 
(fT5|) as a stochastic differential equation 



drj = Xi]dt + ddddW 



(17) 



with the stability parameter A defined in Equation (|15[) . 
The first term on the r.h.s. of Equation (|T7| describes 
the deterministic part of time evolution and the second 
term the stochastic part with dW representing a Wiener 
process with zero mean and variance dt (27| The 
parameter is a measure of the strength of the noise 
term. The solution of the Equation p7|) is given by fl7l . 

M 



V(t) = 77(0)e 



,\(t-s) 



dW(s) 



(18) 



where the temporal variable s is integrated over time to 
time t. In the limit of t — > a, the autocorrelation /^(l) 
and the variance a 2 of ij(t) are determined by the stability 
parameter A [l?], El as 



a 2 = -^ 
" 2IAI 



(19) 



(20) 



with A — ve. The variance a 2 diverges and the lag-1 auto- 
correlation /^(l) reaches its maximum value at the bifur- 
cation point since the stability parameter A is zero at this 
point. Rising variance and autocorrelation thus provide 
early signatures of impending regime shifts. The CSD 
close to the bifurcation point implies that the system's 
intrinsic rates of change are decreased so that the state of 
the system at time t resembles closely the state at time 
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Figure 7. For Model 1, the variation of the stability parame- 
ter A (Equation (|15|l '). the return time Tr, the lag-1 autocor- 
relation p v (l) and the variance (Equations (|19p and (J2S1)) 
as a function of the bifurcation parameter. The plots exhibit 
characteristic features at the bifurcation points k a ,i = 1.344 
and k a< 2 = 5.176. The parameter values are kd = 1, k s — 15, 
k„ = 0.03, Jo = 0.01, Ji = 1 and a d = 0.25. 
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Figure 8. The computed data points in the lnT^ versus (a) 
In \k a — k a ,i\ and (b) In \k a — k a .2\ plots are fitted by straight 
lines (a) y — A + Bx and (b) y = C + Dx respectively. 
The parameter values are the same as in the case of Figure 
□ The values of A and C are A = 3.53401 ± 0.00068 and 
C = 4.35482 ± 0.00049. The exponents have values very close 
to -| (B = -0.49752±0.00013 and D = -0.49901±0.00099). 



t — 1 . This increased memory is measured by the lag-1 
autocorrelation function. Also, because of a flatter basin 
of attraction close the transition point, a given perturba- 
tion brings about a greater shift in the system's state, i.e., 
an increasing variance. 

Two other early signatures of a regime shift which are 
not related to the CSD are skewness and flickering 0, [B| . 
The skewness 7 of a probability distribution P(x) is a 
dimensionless measure of its degree of asymmetry and is 
defined as 



7 



J{x- nfP{x)dx 



(21) 



where \x and a are respectively the mean and standard 
deviation of P(x). The variance a 2 , the coefficient of vari- 
ation C v and the third moment M3 of P{x) are expressed 
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Figure 9. For Model 1, the plot of r = versus the bifurca- 
tion parameter k a where 7\ and T2 are the MFPTs. The ratio 
r is seen to diverge at the lower bifurcation point. Inset shows 
the time-series data of x(t) versus t obtained by solving the 
LE. The straight lines are drawn at the steady state points 
x%, x u , xf obtained by solving the deterministic rate equa- 
tion, Equation ((4]). The parameter values are kd = 1, k s — 15, 
k p = 0.03, Jo = 0.01, Ji = 1 and d 2 = 0.25. 



as 



M, 



C v = 



(x — fi) 2 P(x)dx 
(x - n) 3 P(x)dx 



(22) 



The skewness of a probability distribution increases as 
the bifurcation point is approached. This is because of 
the asymmetry in fluctuations with a system exhibiting 
greater amplitude deviations in the direction of the state 
it is fated to switch to than in the other direction. The 
phenomenon of flickering is observed in the region of bista- 
bility with the system switching back and forth between 
the two attractor states before the bifurcation point is 
reached. We propose a quantitative measure of flickering 
which serves as an early signature of regime shift. In the 
region of bistability, the stability landscape has two min- 
ima corresponding to the two stable steady states. Noise- 
induced transitions take the system from one valley to the 
other. The mean first passage time (MFPT) refers to the 
average first exit time from a valley (2(| [27| ■ Let T\ and 
T2, be the MFPTs for the exits respectively from valley 
1 and valley 2. The times are indicative of the amount of 
flickering present in the system. The MFPT T2 becomes 
zero at the bifurcation point where the steady state 2 
loses stability. The ratio r = ^ measures the asymmetry 
in the exit times and diverges as the bifurcation point, 
at which Ti — ¥ is approached. The quantity r thus 
provides an early signature of regime shift. 



III. RESULTS ON EARLY SIGNATURES 

We first present the results on early signatures in the 
case of Model 1. In the region of bistability, the steady 
state probability distribution P s t(x) is bimodal with two 
distinct peaks. P s t{x) may be obtained via a solution of 



the FP equation (Equation (fT2"j) ) or obtained by a nu- 
merical solution of the LE (Equation©). Figure 0] shows 
an example of a bimodal distribution with two distinct 
peaks at x = x\ and X3. The inset of the figure shows the 
corresponding stochastic potential profile (Equation HU) 
with xi denoting the location of the hill and x\ , X3 denot- 
ing the minima of the left and right valley respectively. 
We find the normalised probability distributions, P™(x) 
and P™(x), corresponding to the low and high expression 
states respectively, using the following procedure: 

1. From the stochastic potential profile 4>f(x), deter- 
mine the position of the hill at xi and compute 
Pst(x2). The point xi corresponds to the minimum 
of the probability distribution between its two peaks 
at .Tiand X3. The point xi is chosen to be the left 
cut-off point, x^, for P™(x). The right cut-off point, 
x R , of P™(x) is obtained as a solution of the equa- 
tion 



P at {x) = P st (x£), x > x 3 



(23) 



2. If P s t(0) < P s t(x2)i then the lower-cutoff point, 
x l L , of P™(x) is determined from the solution of the 
equation 



P st (x) = P st {x l R ), x <Xi 



(24) 



where x l R , the right cut-off point of P™(x), is given 
by x l R = x\. If P st (0) > P st (x 2 ), x l L =Q. 

Other cut-off procedures, e.g., that in Ref. @, have been 
proposed to isolate the low and high expression proba- 
bility distributions but the general results on the early 
signatures are qualitatively similar. 

In the case of Model 1, we first consider only an ad- 
ditive noise in the LE (Equation ((TJ), gi(x) = 0). The 
steady state probability distribution, P jt (x), in the re- 
gion of bistability is obtained from Equation (|12[) putting 
gi(x), gi(x') = 0. Following the prescription already 
given, one determines the cut-off points of the distribu- 
tion Pp(x). The normalised distributions are obtained 
from 



F«P st (x)dx 



for x L < x < x R 



J**P st (x)dx 



for 1? < x < x 



R 



(25) 



(26) 



In Equations ([2~5"j) and (|2l)|) . P st {x) is not normalised. 
With a knowledge of the normalised probability distribu- 
tions P™(x) and P"(x), one can compute the skewness 7, 
variance a 2 , third moment M3 and the coefficient of vari- 
ation C v using Equations (|2"Tj) and (p?2"|) . Figure [5] shows 
plots of these quantities (solid lines) versus the bifurca- 
tion parameter k a for the probability distribution P™(x), 
i.e., considering the system to be in the high expression 
state. The sudden regime shift in the deterministic case 
occurs at the lower bifurcation point k a ,i = 1.334. The 
parameter values used for the computation are: kd = 1, 
k s = 15, k p = 0.03, J = 0.01 and Ji = 1. The strength 
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Figure 10. The phase diagram of Model 2 in the b k — b s plane. 
It shows two regions of monostability, one region of bistability 
and one region of excitability. The boundaries shown as solid 
lines correspond to saddle-node (SN) bifurcations. Lines 1 and 
2 mark the values of the parameter b 3 for which computations 
are carried out. The other parameter values are the same as 
in Ref. [gjj, namely, a k = 0.005, k = 0.2, ki = 0.222, n = 2 
and p = 5. 



d,2 of the additive noise is = 0.25. One finds that all the 
four quantities I7I, <r 2 , M3 and C v increase as the lower 
bifurcation point is approached thus providing early sig- 
natures of a regime shift. The quantities, however, reach 
their maxima before the deterministic bifurcation point is 
reached and then start decreasing. The quantities, though 
providing early signatures, cannot provide knowledge of 
the bifurcation point. We next include an additional mul- 
tiplicative noise term in the LE (Equation ([?])) with 



9i{x) 



k' 



(27) 



The multiplicative noise is associated with the protein 
synthesis rate constant Ji in the active state, i.e., Ji — > 
J 1 +e(t). With both the additive and multiplicative noise 
terms present, the steady state probability distributions 
Pp(x) and Pu(x) are computed following the procedure 
already described. The dotted curves in Figure [3] show 
the variations of I7I, a 2 , M3 and C v , associated with the 
probability distribution P"(.t), as a function of the bifur- 
cation parameter k a . The parameter values are the same 
as before with d\ , the strength of the multiplicative noise 
term, having the value d\ = 0.25. The cross-correlation 
coefficient Xd in Equation © is taken to be zero. Figure 
[5] shows how the skewness of P™ (x) changes as a func- 
tion of the additive noise strength di for the parameter 
values k a = 1.5, kd = 1, k s = 15, k p = 0.03, Jo = 0.01, 
J\ = 1 and d\ = (no multiplicative noise). The value 
k a = 1.5 is greater than the value of = 1.344. The ris- 
ing skewness is thus a signature of noise-induced regime 
shift. In the case of the steady state probability distri- 
bution P"(x), one obtains early signatures of the upper 
bifurcation point k a ^ 2 = 5.176 similar to the ones shown 
in Figure [5] and [6] in the case of P"(x), though the quan- 
titative measures exhibit less prominent variation. 

Considering Model 1, we next calculate the stability 
parameter A (Equation (|15[) h with f(x) given by the ex- 
pression on the r.h.s. of Equation (|4]), in the region of 
bistability. Knowing A, the return time Tr(~ pjy) as well 



as the lag-1 autocorrelation /0^(1) and the variance ct 2 
(Equations (fT9")l ) and (|2H|) ) are also determined. Figure [J] 
shows the variation of A, Tr, Pr)(\) and cr 2 as a function of 
the bifurcation parameter k a and for the parameter val- 
ues kd = 1, k s = 15, k p = 0.03, Jo = 0.01, J\ — 1 and 
Od = 0.25. The stability parameter A becomes zero at the 
lower (k a ,i = 1.344) and the upper (fc Qi2 = 5.176) bifur- 
cation points as expected. The associated return time Tr 
diverges at the bifurcation points, a characteristic feature 
of the CSD. The lag-1 autocorrelation p v (l) reaches the 
maximum value and the variance tr 2 diverges at the bi- 
furcation points. These quantities are good indicators of 
regime shifts and carry distinct signatures of the bifurca- 
tion points. The return time Tr is known to satisfy a gen- 

_ 1 

eral scaling law Tr ~ \B — E>i\ 2 where B and Bi stand 
for the bifurcation parameter and point respectively j3l| . 
In Figure HI the scaling is demonstrated for the bifurca- 
tion parameter k a and the bifurcation points fc a i = 1.344 
and k a ,2 = 5.176. The exponent is very close to — \ in 
each case. 



Figure IHl shows the variation of r = T1/T2 versus the 
bifurcation parameter k a . The ratio is seen to diverge at 
the lower bifurcation point k a> i = 1.344. The parameter 
values are kd = 1, k 8 = 15, k p = 0.03, Jo = 0.01, Ji = 1 
and di = 0.25 (only additive noise is considered). The 
time-series data shown in the inset of Figure [9] is obtained 
via numerical solution of the LE (Equation ([7])) using the 
algorithm described in Ref. [35| . In the limit of large 
times, the steady state is assumed to be reached. The 
solid lines in the inset mark the stable expression states 
Xi and x s h and the dotted line corresponds to the unstable 
steady state x u . The values are obtained from a solution 
of the deterministic rate equation. The MFPTs T\ and T2 
are computed using the method outlined in Ref. [36[ . Let 
us consider a bistable potential with the stable steady 
states at x\ and X3 (x\ < X3) which are separated by 
an unstable steady state at X2, termed the barrier state 
(Figure E]). The MFPT, T(x;a,b) is the average time 
of the first exit from the interval (a, b) and satisfies the 
equation [13, HI 



- 1 = A(x) 



dT(x) 
dx 



2 (X) dx 2 



(28) 



where A(x) and B(x) appear in the associated FP equa- 
tion (Equation ©). The MFPT T(x x )(= Tj) for exit 
from the basin of attraction of the stable steady state at 
x\ is obtained as a solution of Equation (|2"8")l with the 
interval (a, b) = (0, X2) and boundary conditions given by 



T'(a;a,b) = and T(b;a 7 b) = 



(29) 



The prime denotes differentiation with respect to x, with 
reflecting and absorbin g b oundary conditions prevailing 
at a and b respectively |27l . [H| . Following the same pro- 
cedure, the MFPT T(x 3 )(= T 2 ) for exit from the basin 
with the stable steady state at X3 can be calculated from 
Equation (|28[) . The interval now is (a, b) — (2:2,00) with 
X2 and 00 serving as absorbing and reflecting boundary 
points respectively, i.e., 



T(a; a,b) = and T'(b; a,b)=0 



(30) 
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Figure 11. For Model 2, the variation of X max and Tr as a 
function of b s . (a), (b): the computations are carried out along 
the line bt^i = 0.0875. (c), (d): the computations are carried 
out along the line bk,2 = 0.15. The parameter values are the 
same as in the case of Figure ITOT 



We next consider Model 2 the dynamics of which are 
governed by the set of two equations (JSJ) and ((BJ). In Ref. 
[25| , the phase diagram of the model in the b k — b s plane 
has been determined which has four different regions: (i) 
monostability with only one fixed point, (ii) bistability 
with three fixed points, two stable and one unstable, (iii) 
excitability involving three fixed points only one of which, 
corresponding to low ComK value, is stable. The compe- 
tent fixed point (high ComK level) has the characteris- 
tic of an unstable spiral and the mid-ComK fixed point 
is a saddle point, (iv) only one fixed point exists which 
is unstable. The system exhibits limit cycle oscillations 
between the mid-ComK and high-ComK levels. For the 
purpose of our study, a part of the bk — b s phase diagram 
has been recomputed and shown in Figure [TDJ The dia- 
gram shows three different regions: monostable, bistable 
and excitable. The monostable region is again of two 
types: monostable low (low ComK level) and monostable 
high (high ComK level). The boundaries between the 
regions depicted by solid lines correspond to the saddle- 
node (SN) bifurcation Two vertical lines 1 and 2 are 
drawn in the phase diagram at the points bk,i = 0.0875 
and bk,2 = 0.15 respectively. Line 1 intersects the phase 
boundaries at the three points b 8> i — 0.7799, b Sj2 = 0.7868 
and b s .3 — 0.8094. Line 2 has two points of intersection: 
b Sll = 0.6175 and b s , 2 = 0.7504. Figure HU shows the 
plots of Xmax, the real part of the dominant eigenvalue 
of the Jacobian matrix computed from Equations ([5J and 
©, and the return time Tr = — r versus b s along line 
1 ((a) and (b)) and along line 2 ((c) and (d)). For a spe- 
cific stead y st ate, the Jacobian matrix J has the following 
structure |13j : 



J 



fll fl2 
721 1 22 



(31) 



ss 



The suffix SS stands for steady state, i.e., the matrix 
elements are to be computed at the fixed point (x*,y*). 



A matrix element = (i = 1, 2, j = 1, 2) where 

bk X! 



fi(xx,x 2 ) = a k + — 



fc£ + x% 1 + x x + x 2 



f2(Xl,X 2 ) 



x 2 



l + ixi/klY l+X!+X 2 



(32) 



(33) 



with x\ = K, x 2 = S, n = 2 and p = 5. The param- 
eter values are the same as in Ref. [25j : a k = 0.005, 
fco = 0.2, and k\ = 0.222. The eigenvalue of J are Ai and 
A2 and X m ax is the real part of the dominant eigenvalue. 
In Figure II H one notes that X max becomes zero at the 
SN bifurcation points as expected and the correspond- 
ing return time Tr diverges at a bifurcation point. The 
shaded regions in the Figure denote the regions of bista- 
bility in which the two stable steady states correspond to 
low and high ComK levels respectively. The level of ComS 
is anticorrelated with that of ComK. The solid (dotted) 
lines in Figure [11] are associated with the stable steady 
states representing low (high) ComK levels. Along line 2 
and at the bifurcation point 6 Sj i = 0.6176 (Figures ITTT c') 
and (d)) the stable steady state corresponding to the high 
ComK level loses stability (X max = 0, Tr diverges). At 
the boundary point b S:2 = 0.7504, the state representing 
the low ComK level loses stability. A steady state is sta- 
ble if the real parts of both Ai and A2 are negative. At the 
point 6 s .i = 0.7799 along line 1 (Figures [TTJa) and (b)), 
there is no loss of stability of a steady state and hence no 
CSD with diverging Tr is observed at this point. The low 
ComK level of the monostable region continues to remain 
stable as one enters the region of excitability at the point 
6 Sj i along line 1. At the point b St2 = 0.7868, the system 
enters the region of bistability from a region of monosta- 
bility. The high ComK stable steady state loses its sta- 
bility at this point with X max = and a divergent Tr 
(dotted branch in Figures ITTT a) and (b)). The low ComK 
state loses stability at the point b s ,3 — 0.8094 when the 
system passes from a region of bistability to a region of 
monostable high ComK level. The eigenvalue X max = 
and Tr diverges at this point. The main point to note 
from the results is that there is no CSD and diverging 
return time Tr at the transition point b s _\ along line 1 
between the regions of monostability and excitability. On 
the other hand, at the point b s _i along line 2, separat- 
ing the regions of monostability and bistability, one of 
the expression states, namely, the high ComK state, loses 
stability with a divergent Tr. 



IV. CONCLUDING REMARKS 

In this paper, we investigate the early signatures of 
sudden regime shifts occurring at the bifurcation points 
associated with a fold-bifurcation model. The basic con- 
cepts and methodology are applicable to a general class of 
bifurcation phenomena occurring in diverse systems. Our 
focus, however, is on obtaining quantitative estimates of 
the early signatures of regime shifts in the gene expres- 
sion dynamics of competence development in B. subtilis. 
In the case of Model 1, which is a one variable model, 
the formalisms based on the LE, the FP equation and a 
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stochastic differential equation (Equation (fTT)) ) are uti- 
lized to compute the different quantities. The early indi- 
cators of regime shifts, which include the CSD, variance, 
lag-1 autocorrelation, skewness and the ratio r of MFPTs, 
arc shown to exhibit distinctive variations as a bifurcation 
point is approached. Though the literature on the early 
signatures of regime shifts in ecosystems, financial mar- 
kets and complex diseases is extensive [l|-|8[ , the issue has 
not been systematically addressed in the case of gene ex- 
pression dynamics, a fundamental activity in the living 
cell. Our study is the first in this direction and illustrates 
how the early signatures provide knowledge in advance of 
an impending regime shift. Some of these signatures may 
also provide clues on the physical principles underlying 
gene expression dynamics. As already pointed out, X max 
and Tr do not exhibit any distinctive features when the 
system switches between a region of monostability (low 
ComK level) and a region of excitability (Figure [TlTa)). 
This is in contrast to the case when the system passes 
from a region of bistability to one of monostability. The 
quantity X ma xj associated with the steady state which 
loses stability, becomes zero at the bifurcation point and 
the return time Tr diverges. Experiments detecting the 
CSD, as the bifurcation point is approached, are difficult 
to carry out. Observation of the CSD in recent experi- 
ments [32l . l33l | provides pointers for carrying out further 
such experiments. The CSD experiment, if designed in 
the appropriate manner, would be able to distinguish be- 
tween the bistability versus excitability paradigm in the 
case of competence development in B. subtilis. In the case 
of Model 2, we have not obtained quantitative estimates 
of the other signatures like rising variance and lag-1 auto- 
correlation function, skewness and flickering because the 
FP and stochastic differential approaches are computa- 
tionally challenging for more than one variable. Also, the 
rate equations (Equations (J5J) and (JB])) provide an effec- 
tive description of the gene expression dynamics based 
on various approximations and simplifications. A more 
faithful representation of the dynamics would be in terms 
of a larger number of rate equations involving more than 
two variables [25|. We defer the computations on early 
signatures in the case of Model 2 to a future study. 



One utility of obtaining knowledge of the early signa- 
tures is that specific risk aversion strategies may be devel- 
oped if the sudden transition is to a regime which is un- 
desirable due to various considerations. Examples of such 
regime shifts include asthma attacks Q , epileptic seizures 
[H and the sudden deterioration of complex diseases [l8| . 
One may add the development of persistence of pathogens 
like M. tuberculosis in the human lung granulomas to the 
list. Recent experiments H3, H3 on a sister species M. 
smegmatis provide evidence that a fraction of the my- 
cobacterial population (the population of persisters) is 
able to survive under nutrient depletion. To achieve this, 
the mycobacteria adopt the strategy of generating phc- 
notypic heterogeneity in the form of two distinct subpop- 
ulations. The concentration of a key regulatory protein 
Rel is high in one subpopulation (which develops into the 
persister subpopulation) and low in the other. In the 
subpopulation with high Rel level, the stringent response 
pathway is initiated which help the mycobacteria to avoid 
death and adapt to nutrient depletion. The experiments 
(37l . |38| show that the principle underlying the develop- 
ment of phenotypic heterogeneity is based on bistabil- 
ity and noise-induced transitions between the expression 
states corresponding to low and high Rel levels. The prob- 
lem of persisters is that they are not killed by antibiotic 
drugs and wait for the opportune moment to restart an 
infection [39| . Early signatures of a regime shift from the 
normal to the persistent state would help in developing 
measures preventing the switch to persistence. Similar 
studies could be carried out on other systems in which bi- 
furcation phenomena are responsible for the development 
of heterogeneity, choice of cell fate [l9| or regime shifts 
leading to new types of dynamical behaviour. 
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